План митапа

  • Типы пространственных данных;

  • Особенности геоданных;

  • R и пространственный анализ;

  • Практика - находим геоданные разных форматов и учимся с ними работать;

  • Что можно почитать и посмотреть по теме.

 

1. Типы пространственных данных

 

а. Векторные данные

  • Геометрии состоят из точек. Конвенционально x = широта, y = долгота. Кроме этого могут использоваться z = высота, M = мера изменчивости точки;

  • Виды геометрий:

+++ GEOMETRYCOLLECTION +++

  • Форматы записи данных:

    Well-Known Text (WKT) - Запись точек вектором из геометрий

    Well-Known Binary (WKB) - Запись координат в бинарных значениях. Используется в базах данных (т.к. увеличивает скорость работы с данными), но нечитаем и непонятен для человека.

Посмотрим в первом приближении:

library(sf)

library(rnaturalearth)
world <- ne_countries(returnclass = "sf")
world <- world[1:5,c(4, ncol(world))]

# Для иллюстрации
world_1 <- world
world_1$geometry <- as.character(world_1$geometry)
as.data.table(world_1)

 

b. Растровые данные

  • Матрица значений пикселей пространственной области. В растрах хронят космические снимки и базы различных геологических данных. Также как и вектор - предназначены для определённого круга задач;

  • Основной формат: .tif

  • Основные библиотеки для обработки в R: raster и stars

Сравним:

На этом и последующих занятиях мы будем говорить только о векторных данных.

 

Главное, нужно помнить, что векторные пространственные данные - это (чаще всего) просто точки на плоскости! Создадим собственные пространственные данные:

p1 = st_point(c(7,52)) # Рисуем точку 1
p2 = st_point(c(-30,20)) # Рисуем точку 2
sp_obj = st_sfc(p1, p2, crs = 4326) # Преобразуем в пространственные данные
plot(sp_obj)

# Вы восхитительны!

 

2. Особенности геоданных

Вернёмся к векторным геометриям:

world

 

Coordinate Reference Systems (CRS) / Пространственная привязка и её проекции

Что это такое?

Для иллюстрации обратимся к прекрасной работе Тасс.

 

Картинки чтобы окончательно закрепить идею проекции.

 

 

 

Как посмотреть проекцию в CRS?

st_crs(world)
## Coordinate Reference System:
##   User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

 

 

На самом деле всё не так сложно.

Пространственный объект можно перевести в другую CRS сразу и быстро, задав через st_transform() нужную проекцию одним из четырёх форматов (просто гуглите то, что нужно для вашего пространства):

  • Код EPSG (например, “4326”)
  • Формула PROJ4 (например, “+proj=longlat +datum=WGS84 +no_defs”)
  • Строчка WKT
  • Просто переносом проекции с другого пространственного объекта ( df1 <- st_transform(df1, st_crs(df2)) )

Например:

Afg = world[1,]
plot(Afg)

# Change crs
Afg <- st_transform(Afg, 2264)
plot(Afg)

# Change crs
Afg <- st_transform(Afg, "+proj=longlat +datum=WGS84 +no_defs")
plot(Afg)

В геоданных есть ещё множество подводных камней, но этот - ключевой и его необходимо запомнить.

 

3. R и пространственный анализ

Почему R?

 

 

Библиотеки для работы с геоданными в R

(несть им числа, но основные)

  • sf - state-of-the-art пространственного анализа в R.

Немного истории:

 

4. Практика

Что мы будем делать - делать карту мороженного в Израиле! :)

Я исхожу из того, что вы знаете dplyr или data.table на базовом уровне, но если у вас вызывает затруднения именно манипуляции с данными - обязательно задавайте вопросы!

Мы отработаем работу с данными на трёх “классических” форматах для векторных данных: csv, geojson, shape-file.

 

A. Датасет с мороженками в csv

Загрузим нужные библиотеки:

library(sf)
library(geojsonsf)
library(geojsonio)
library(readr)
library(data.table)

Загрузим данные. К сожалению, в R есть известная проблема с кодировками и чтобы долго не решать проблему с крокозябрами - уберём колонки названий точек мороженного на иврите:

isr_icecream <- fread('https://raw.githubusercontent.com/AmitLevinson/Datasets/master/golda/golda_locations.csv', encoding = 'UTF-8')
isr_icecream[,1:2 := NULL]
isr_icecream

Превратим в sf тип данных:

isr_icecream <- st_as_sf(isr_icecream, dim = "XY", remove = T, na.fail = F, coords = c("lon", "lat"), crs = "+proj=longlat +datum=WGS84 +no_defs")
isr_icecream
plot(isr_icecream)

 

B. Границы Израиля в Shape-file

Идём сюда, ищем Израиль и загружаем архив. Что мы видим внутри?

.shp – главноеый файл с геометриями

.dbf

.shx

.prj

(но все они важны!)

По этой причине, когда вам нужно переслать shape-file отправляйте архив из всех четырёх файлов.

Загрузим его и отобразим:

isr_border <- st_read('/cloud/project/ISR_adm/ISR_adm0.shp')
## Reading layer `ISR_adm0' from data source `/cloud/project/ISR_adm/ISR_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 34.26801 ymin: 29.49708 xmax: 35.90094 ymax: 33.36403
## CRS:            4326
isr_border
plot(st_geometry(isr_border))

 

C. Точки городов Израиля в GeoJson

Json с геоданными. Достаточно тяжёлый с точки зрения размещения данных, но часто встречающийся формат.

И здесь стоит отвлечься и ознакомиться с “Википедией”" от мира пространственных данных - Open Street Map.

Если это “Википедия”, то с неё можно выкачать данные? Да. По тэгам…

  …  

Тэги городов.

Скачаем GeoJson с городами. Загрузим его. Внутри будет огромное количество колонок, которые сейчас нам не нужны. Оставим первую и последнюю:

isr_city <- st_read('/cloud/project/city.geojson')
## Reading layer `city' from data source `/cloud/project/city.geojson' using driver `GeoJSON'
## Simple feature collection with 51 features and 225 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 33.80459 ymin: 30.61197 xmax: 36.56875 ymax: 32.92817
## CRS:            4326
isr_city <-  isr_city[,c(1,ncol(isr_city))] 
isr_city
plot(isr_city)

 

Альтернативный способ через пакет osmdata (спасибо Philipp Pyshny за наводку)

library(osmdata)

admin_osm = opq(bbox = 'Saint-Petersburg') %>%
  add_osm_feature(key = "admin_level", value = '5') %>%
  osmdata_sf()

# Но есть проблемы с кодировкой
iconv(admin_osm$nodes$tags, from="UTF-8", to="UTF-8") 
iconv(admin_osm$nodes$tags, from="UTF-8", to="cp1251")
admin_osm[["osm_multipolygons"]]

 

Объединяем слои и строим свою первую карту!

 

# Устанавливаем одинаковые проекции для всех объектов: 
isr_border <- st_transform(isr_border, 4326)
isr_city = st_transform(isr_city, st_crs(isr_border))
isr_icecream = st_transform(isr_icecream, st_crs(isr_border))

# Делаем карту (st_geometry() - вывести чистую геометрию, без других колонок или значений в объекте)
{plot(isr_border %>% st_geometry())
plot(isr_city, col = 'red', add = T)
plot(isr_icecream, col = 'blue', add = T)}

 

Ура :)

 

Теперь сохраним данные (здесь есть свои подводные камни):

# write_sf(isr_border, "isr_border.shp", append = T, layer_options = "ENCODING=UTF-8")

Загрузим в качестве проверки и убедимся, что всё работает

# double_isr_border <- st_read("isr_border.shp")
# plot(st_geometry(double_isr_border))

 

Бонус

Построим интерактивный график с помощью очень простой по своему синтаксису библиотеки mapview. Туториал по ней можно найти здесь. mapview предназначен для карт небольшого размера.

ЗАПУСКАТЬ ОСТОРОЖНО!!!

(на Windows возможны вылеты RStudio)

library(mapview)
mapview(isr_border) + 
  mapview(isr_city %>% st_geometry(), col.regions = 'red', legend = FALSE) +
  mapview(isr_icecream %>% st_geometry(), col.regions = 'blue', legend = FALSE)

 

5. А где брать данные?

Open street map

naturalearthdata.com/downloads/

download.geofabrik.de/

#30DaysMapChallenge

Один из многочисленных обзоров челленджа

и много, много других ресурсов…

 

6. Что почитать?

Pebesma E., Bivand R. Spatial Data Science. 2020.

Lovelace R., Nowosad J., Muenchow J. Geocomputation with R. 2021

Dorman M. Using R for Spatial Data Analysis. 2021.

Простое гугление даёт как минимум 10 книг о пространственном анализе в R.